An unconstrained approach to systematic structural and energetic screening of materials interfaces

From grain boundaries and heterojunctions to manipulating 2D materials, solid-solid interfaces play a key role in many technological applications. Understanding and predicting properties of these complex systems present an ongoing and increasingly important challenge. Over the last few decades computer simulation of interfaces has become vastly more powerful and sophisticated. However, theoretical interface screening remains based on largely heuristic methods and is strongly biased to systems that are amenable to modelling within constrained periodic cell approaches. Here we present an unconstrained and generally applicable non-periodic screening approach for systematic exploration of material’s interfaces based on extracting and aligning disks from periodic reference slabs. Our disk interface method directly and accurately describes how interface structure and energetic stability depends on arbitrary relative displacements and twist angles of two interacting surfaces. The resultant detailed and comprehensive energetic stability maps provide a global perspective for understanding and designing interfaces. We confirm the power and utility of our method with respect to the catalytically important TiO2 anatase (101)/(001) and TiO2 anatase (101)/rutile (110) interfaces.

In Supplementary figure 3 we compare a displacement scan at α = 0° calculated when using: i) our adopted relaxation protocol (i.e. atoms relaxed within a radius of 1.7 nm and thickness of 1 nm with all outer atoms fixed to bulk positions), and ii) disks with all atoms fixed to bulk position. Clearly, allowing inner interfacial region to relax slightly lowers the calculated energy in a systematic way.
To show this effect more clearly, we take three characteristic points from these maps, corresponding to: an Esep maximum (max), Esep minimum (min) and an intermediate Esep between these values (between).
We then follow the change in energy of each with respect to the percentage of the relaxed interfacial area (i.e. with respect to the maximum corresponding to the full 2.5 nm radius). We can see in each case that Esep systematically decreases in a similar linear fashion with respect to the percentage of relaxed interfacial area employed. This implies that all relative interface energetic stabilities within a displacement scan are relatively unaffected. From this tendency we can also estimate that the Esep values that we obtain using a relaxed region of 1.7 nm radius are systematically overestimated by 0.2 J/m 2 relative to a fully relaxed system. Atom key -Ti -blue, O -red.

Supplementary note 6
For any selected displacement we can perform a systematic twist angle scan. In the main text we use an angular increment of 5° for these scans. Below we show that the results are not significantly changed if we increase the resolution of the scan by five times by using a 1° increment indicating that 5° is sufficient for an initial twist angle scan.

Supplementary note 7
The calculated a and c lattice parameters obtained from an IP-based optimization of the anatase TiO2 bulk crystal structure were found to be 3.770 and 9.570 Å respectively. The corresponding values predicted by periodic DFT calculations with the HSE06 functional were 3.775 and 9.478 Å, respectively. The IP-based and DFT based values are in good agreement with a deviation <1 % for both a and c. Both sets of values are also in good agreement with previously reported DFT data and experimental measurements, a = 3.785 and c = 9.514 Å. An IP versus DFT comparison of the structures of the optimised periodic (001) and (101)  The calculated lattice parameters of rutile (110) obtaining with our employed IPs are a = 3.013 Å and b = 6.337 Å, which are again in good agreement with the corresponding DFT results, 3.01 Å and 6.56 Å respectively.

Supplementary note 9
The coincidence site lattice (CSL) approach is traditionally applied to grain boundaries (GBs) of densely packed monoelemental systems (typically metals). For a pure twist GB one can first align the lattice parameters of the two grains and then twist one lattice with respect to the other until one finds specific angles with periodic lattice site coincidences. Although CSL does not consider atoms, the idea is based on finding GBs with a good atomic fit and are which are thus expected to have higher energetic stabilities. CSL has often proven to be a useful geometric guide to low energy metal GBs.
In our work, we directly evaluate the relative energetic stabilities of interfaces between dissimilar surfaces of a semi-ionic compound in a non-periodic system (i.e. not a typical GB system). Simple geometric measures such as CSL for understanding interface stability have been reported as being lacking with respect to more direct atomistic models (e.g. "It is concluded that no general and useful criterion for low energy can be enshrined in a simple geometric framework. Any understanding of the variations of interfacial energy must take account of the atomic structure and the details of the bonding at the interface." -A. P. Sutton, R. W. Balluffi (1987), Overview no. 61: On geometric criteria for low interfacial energy, Acta Metallurgica, 35, 2177). In our study we wanted to use a CSL-type measure that was well tailored for our systems and thus less 9 susceptible to such criticism. In particular, the relatively high complexity of our systems means that a traditional CSL approach has limited applicability: 1.) When considering heterojunctions, traditional CSL does not give a unique prescription for finding the relative lattice displacement from where to start twisting. Therefore, we first perform a scan of rigid body displacements to obtain a low energy starting point for our twist angle scans. For consistency, we use these displacements for both our CS scans and our energy evaluations with respect to twist angle. Note that these displacements are determined independently of the underlying lattice sites for the respective surfaces.
2.) Traditional CSL is purely geometric and does not consider atoms nor discriminate between atom types/charges. Thus, CSLs for atomically heterogeneous systems can yield highly repulsive and energetically unstable interfaces. To take into account such cases, we give each atomic site a sign and a local spatial region (Gaussian distribution) and calculate the overlap between cation/anion sites in our coincidence measure.
This means that we go beyond a purely lattice based approach to a measure which uses atoms positions when evaluating CSs. We note that CS approaches to rationalise interfaces using atoms positions are not new and were first proposed and used effectively in: Kronberg, M. L. and F. H. Wilson (1949), "Secondary recrystallization in copper", Trans. Met. Soc. AIME, 185, 501.
Considering the above, to provide an atomic CS-based evaluation of our interfaces we assigned to each atom (N) having coordinates (x0,y0) at the interface between the two disks, an isotropic gaussian function of the type: The Gaussian width (s) was set to be 0.1, which was found provide stable smoothing varying results. In this way, two atoms belonging to two different disks, but having the same x and y atomic coordinates at the crystallographic plane of the interface, would have the maximum overlap, calculated as: We note that within this approach the positioning along the non-periodic direction z is not considered. In other words, this method provides most likely overlap atomic configurations in the plane defining the interface, and the z-coordinates will be given by the ideal atomic bond distances between the atoms belonging to the two disks.
In addition to simple atomic overlap, to account for the ionic nature of TiO2 we use a signed overlap measure to account for cation-anion and cation-cation or anion-anion matching by assigning a positive or negative sign to atoms according to the sign of their formal charge. In such a scheme, a perfect overlap (coincidence) between two cations or two anions gives rise to a value of 1. Conversely, the perfect overlap (coincidence) between a cation and an anion, results in an overlap equal to -1. Therefore, each possible overlap range can span from -1 to +1. For each applied displacement and rotation in our double disk scans we calculate the signed overlap integral for all atoms at the interface as:
Finally, we normalize the overall integral by the number of interface atoms to provide a final value between -1 and +1. The resulting signed atom-based CS maps are compared with the relative energy maps with respect to scanning a range of relative displacement and angles in Supplementary figure 10.